arXiv:1509.01696v5 [math.DS] 23 Sep 2016 


Early-warning indicators for rate-induced tipping 
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A dynamical system is said to undergo rate-induced tipping when it fails to track its quasi-equilibrium state 
due to an above-critical-rate change of system parameters. We study a prototypical model for rate-induced 
tipping, the saddle-node normal form subject to time-varying equilibrium drift and noise. We find that both 
most commonly used early-warning indicators, increase in variance and increase in autocorrelation, occur 
not when the equilibrium drift is fastest but with a delay. We explain this delay by demonstrating that the 
most likely trajectory for tipping also crosses the tipping threshold with a delay and therefore the tipping 
itself is delayed. We find solutions of the variational problem determining the most likely tipping path using 
numerical continuation techniques. The result is a systematic study of the most likely tipping time in the 
plane of two parameters, distance from tipping threshold and noise intensity. 
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The notion of tipping describes the phenomenon 
that at certain critical levels or rates the output 
of a system changes disproportionately compared 
to the change in input. Two particular exam¬ 
ples in climate science are the possible collapse of 
the Atlantic Meridional Overturning Circulation 
(AMOC) due to increasing freshwater input or 
the sudden release of carbon in peatlands due to 
an external temperature increase above a critical 
rate (the compost bomb instabilit. There is an 
ongoing debate, for example, in climate scienc^ 
and ecologjEl whether it is possible to find early- 
warning indicators robustly in time series of sys¬ 
tem outputs before tipping occurs. 

This paper focuses on the case of rate-induced 
tipping, which describes the scenario where a sys¬ 
tem fails to track its equilibrium due to a rapid 
change in parameters. We show that two popu¬ 
lar candidates for early-warning indicators in time 
series, an increase in autocorrelation and an in¬ 
crease in variance, appear to give a delayed warn¬ 
ing signal for tipping. We study the phenomenon 
by looking at the interaction of rate-induced tip¬ 
ping and noise. We find that the most likely time 
for the noise to kick the system over the threshold 
(which is a curve in phase space) is after the point 
where the threshold is closest to the equilibrium. 
We investigate this tipping (and, correspondingly, 
early-warning) delay systematically depending on 
two parameters: the distance of the parameter 
drift speed from its critical value and the noise 
intensity. We find that the delay is larger for 
smaller noise-intensity. 
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I. INTRODUCTION 

Tipping events are often described as sudden, dis¬ 
proportionate changes in output levels caused by small 
changes to input levelJ^. These can be irreversible 
events that have huge, unwanted consequences. There¬ 
fore, the study of early-warning indicators is of great 
interest and so recent research has developed and an¬ 
alyzed early-warning indicators, see Lentonl^ for a re¬ 
view up to 2011 and Williamson and Lentonl^ for refer¬ 
ences to later results. A few currently debated examples 
of complex systems deemed vulnerable to tipping from 
climate science, ecology and financial markets are: the 
abrupt reductions in Arctic summer sea ic^, the col¬ 
lapse of the Atlantic Meridional Overturning Circulation 
(AMOC)l^, the dieback of the Amazon rainforesfP, Aus¬ 
tralian ecosystemJISl^ the light-driven regime shifts in po¬ 
lar ecosystemJ^, the collapse of coral reefs due to global 
warming and ocean acidificatiorP^, and crashes and re¬ 
bounds of financial market^. See also Lenton et al^ 
for a list of policy-relevant tipping elements in the climate 
system. Ashwin et alJ^ identified a few mathematical 
mechanisms behind the observed phenomena, attempt¬ 
ing a classification: 

• Bifurcation-induced tipping (Slow passage through 
a bifurcation) 

• Noise-induced tipping (Transition between attrac¬ 
tors due to random fluctuations) 

• Rate-induced tipping (Failure to track a continu¬ 
ously changing quasi-steady state). 

This paper studies how a system that is close to a 
rate-induced tipping event behaves under the influence 
of additive noise. We look at a prototypical system, 
the saddle-node normal form with additive noise and 
a ramped shift of the equilibrium as proposed by Ash¬ 
win et al^. (A more general definition and further 
properties of rate-induced tipping are given by Ashwin, 
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Perryman, and WieczorekSl.) Two early-warning indi¬ 
cators that are commonly used for bifurcation-induced 
tipping with noise are an increase of autocorrelation and 
an increase in variance in observed time series of system 
outputThe most common argument, why a generic 
output time series of a system approaching bifurcation- 
induced tipping should show an increase in autocorre¬ 
lation and variance assumes that the bifurcation of the 
deterministic part is a saddle-node bifurcation. Far away 
from the bifurcation one can think of the state of the 
system as the position of an overdamped particle at the 
bottom of a slowly softening potential well. Any small 
perturbation or disturbance will relax back to the equilib¬ 
rium with a large decay rat^^^. As the bifurcation is ap¬ 
proached, the potential well will become shallower, that 
is, the decay rate will decreas^^^. Thus, any disturbance 
or perturbation will have an increased and more long- 
lasting effect such that, in the presence of noise, autocor¬ 
relation and variance in observed time series increases. 

In practice, the early-warning indicators are used 
on observational time series data of systems where 
quantitatively accurate models are unavailable, such as 
palaeoclimate temperature and CO 2 proxie^ and lake 
eutrophicatiorP^. In the cited cases the early-warning 
indicators were not used for prediction (as they were 
about events in the past), but as evidence for (or against) 
the presence of und erly ing tipping mechanisms. For ex¬ 
ample, Dakos et used an increase in autocorrela¬ 

tion in a sequence of palaeoclimate time series as evi¬ 
dence for bifurcation-induced tipping while Ditlevsen and 
Johnsen^^ used the absence of the increase in variance 
(and the inconclusive behavior of the autocorrelation) as 
evidence that the Dansgaard-Oeschger events are a case 
of noise-induced tipping. Similarly, the presence of early- 
warning indicators in simulation data from a global cir¬ 
culation model showing a collapse of the Atlantic Merid¬ 
ional Overturning circulation (AMOC) through freshwa¬ 
ter input was used as evidence for a bifurcation-induced 
tipping evenfP. All of the cited studies base their ar¬ 
guments on the knowledge that, close to a bifurcation- 
or noise-induced tipping point the system (after de¬ 
trending) behaves like an Ornstein-Uhlenbeck (OU) pro¬ 
cess. For the OU process one can infer from observed 
autocorrelation and variance the underlying linear decay 
rate, thus, permitting conclusions about the approach (or 
lack of it) of the equilibrium to a saddle-node bifurcation. 
This paper studies the effect of noise on the third mech¬ 
anism from the list by Ashwin et rate-induced tip¬ 

ping, with the goal to aid identification of this type from 
time series. 

In contrast to bifurcation-induced tipping, rate- 
induced tipping is failure of the system to track the 
continuously changing quasi-steady stat^^. Unlike 
bifurcation-induced tipping, at each moment in time 
there exists a stable (quasi-)equilibrium but the rate at 
which this steady state shifts determines whether the sys¬ 
tem tips or not. 


The effect of rate-induced tipping has been described 
only relatively recently. In particular, within climate sci¬ 
ence Wieczorek et a/.® considered a model for carbon 
storage and release in peatland soil, which showed the 
compost bomb instability. In their model an increase in 
temperature above a critical rate results in a release of 
carbon into the atmosphere from combustion of compost 
heaps. A higher CO 2 concentration in the atmosphere, 
creates further warming and thus triggering a positive 
feedback loop within the systenP^. This is an example 
of rate-induced tipping as for every fixed atmospheric 
temperature there exists a globally stable steady state 
but the rapidity of the temperature increase causes sharp 
peaks of carbon release. Other examples of rate-induced 
tipping include the switching off of the AMOC due to 
the rate of increase of CO 2 in the atmospher^^. Scheffer 
et find in a plant-herbivore model the critical rates 
of plant growth causing a rate-induced transition from a 
herbivore controlled state to a vegetated state. 

Rate-induced tipping is not associated to a loss of sta¬ 
bility of equilibrium and thus cannot be explained us¬ 
ing stability theory for equilibrigP^. An appropriate ana¬ 
logue to the “overdamped particle in a softening well” il¬ 
lustration for bifurcation-induced tipping is to think of 
an overdamped particle in a moving well. In contrast to 
bifurcation-induced tipping, the shape of the potential 
well remains constant but instead shifts at varying rates. 
The faster the shift the further the particle drifts away 
from the bottom of the well, up the side and thus closer 
to the saddle and escaping. Hence, there is no change in 
stability of the potential well, only the location of where 
the state is in terms of the potential. As a consequence, 
Ashwin et remarked there is no reason to assume 

why the early-warning indicators: autocorrelation and 
variance can still give useful predictions. 

This paper builds on the work of Ashwin et 
which introduced a prototypical model for deterministic 
rate-induced tipping. We will consider the effect of addi¬ 
tive white noise on this prototype model of rate-induced 
tipping. This models fluctuations/uncertainties that ex¬ 
ist in various systems, for example the climate system. It 
also permits us to study early-warning indicators. The 
aim of this paper is to demonstrate that autocorrelation 
and variance will show an increase. However, this in¬ 
crease occurs with a delay, which is related to a delay in 
the actual tipping. 

The paper is structured as follows: Section |TI| describes 
the basic properties of the deterministic prototype model 
for rate-induced tipping introduced by Ashwin et a/.^. 
Section |HI| explores the apparent delay of the early- 
warning indicators for noise and rate-induced tipping. 
In Section [IV| we set up a boundary-value problem for 
most likely tipping paths, the sequence of continuation 
steps to solve this boundary-value problem are presented 
in Section [V| In Section |Vl] the most likely tipping path 
is discussed for a fixed set of system parameters, and in 
Section [VH| analysis of most likely paths for all relevant 
system parameters using numerical continuation is cov- 
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ered. SectiQn |VIII| discusses results of delay in the context 
of autonomous systems before, Section [represents some 
concluding remarks. 


II. THE DETERMINISTIC BACKBONE — A 
PROTOTYPE FOR RATE-INDUCED TIPPING 


ODE in the (x, A) phase plane (as done by Ashwin 
et alJ^: 

X= fix,X{t)) = {x +X{t)f-1 (4) 

X— ^Xi^XyyiaX '^) ( 5 ) 

Notice that (§ is coupled to 0 , but there is no coupling 
in the other direction. Figure^ displays all qualitatively 


A prototype model for rate-induced tipping was intro¬ 
duced by Ashwin et alJ^. The model is a scalar ordinary 
differential equation (ODE) for the variable x{t) G M: 

X = f{x,X) = {x + xf -1. (1) 

which, is the normal form for the saddle-node bifurcation. 
We have set the normal form parameter equal to 1 w.l.o.g. 
(corresponding to a choice of scale for x and time). The 
ODE 0 has two A-dependent families of equilibria, one 
stable at Xeq (A) = —A — 1 and one unstable at Xeq^(A) = 
—A + 1. The equilibria are separated by a distance of 
2. These families of equilibria Xeq (A) and Xeq^A) form 
straight lines in the (A, x) - plane and will be referred to 
as Wo and Wq respectively (see Figure [^. Equation Q 
is the saddle-node normal form shifted by A, for which we 
assume dependence on time in the form of a ramp (see 





t A 

(a) (b) 



A(t) 


^max 



+ 1 


( 2 ) 


where Amax (distance) and e (speed) are the shape pa¬ 
rameters of the ramp-like shift. The time-derivative of 



FIG. 1: Time profile of shift parameter A(t), equation 
Q, where the black dashed lines indicate the transition 
period for Amax = 3, e = 1 


(c) 


(d) 



(e) 


(f) 


FIG. 2: Time profiles (^,(c),(e) and phase planes 

(b) ,(d),(f) of system 0-(|^ for e < ec - (a),(b), e = ec - 

(c) ,(d) and e> Cc - (eX(f)- Black dashed curves are the 
stable Wq and unstable Wq equilibria in the limit 

e = 0, blue and red curves are the unstable and stable 
manifolds, W'^{S-) and W^(I/+), respectively 
(-^max 3). 


A(t) is of most interest here as this determines the rate 
of shift for the (now quasi-) equilibria Xeq and x^^. The 
time derivative of A is 


dA 

dt 




sech 


Ama,x^^ 


eA(A„ 


A). (3) 


This time derivative reaches its maximum at t = 0 and 
so, for a fixed ramp height Amax, e is directly proportional 
to the maximal rate of shift at t = 0. 

We note that § is also an ODE for A such that the 
prototype model can be considered as a two-dimensional 


different phase portraits possible for 0-0 in panels 
|2d|and[^ The system has 4 equilibria, S- (a saddle) and 
U- (a source) on the A = 0 line, and 5'+ (a sink) and 
(a saddle) on the A = Amax line. The upper and lower 
black dashed lines represent the family of unstable Wq 
and stable Wq quasi-equilibria in the limit e = 0, respec¬ 
tively. The blue curve is the unstable manifold W'^{S-) 
of the saddle S-^ and the red curve is the stable mani¬ 
fold W^(I/+) of the saddle I/+. The panelsand 
show the time profiles for x on the invariant manifolds 
W'^{S-) and W^(I/+) (using the same color coding). 
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One can see that the time profile and the phase portrait 
of the unstable manifold W'^{S -) and the stable manifold 
(Figures and [21^ deviate increasingly from 
the quasi-equilibrium families Wq and Wq for increasing 
e. For small e, is close to Wq , but, for increasing 

e, moves further apart from Wq. The unstable 

manifold W'^{S-) converges for t ^ oc to the stable 
node 5'+ for e < €c. The red curve is the stable manifold 
([/+), which forms a separatrix partitioning the plane 
into two regions. In the region below the separatrix all 
trajectories are attracted to the stable node 5'+, but in 
the region above the repelling stable manifold 
all trajectories will escape to Too in finite time. Notice 
that the two manifolds and are closest 

at A = Amax/2, when the time-derivative A (equation 
is at its maximum. This is due to the reflection symmetry 
within the system 0, (§ 


X — Xn 

A-Ae 


Xn — X 

Ac-A 


around the point (xc, Ac) = (—1.5,1.5). 

At a critical e, denoted ec, and W^{S-) form 

a heteroclinic connection between the two saddles S- and 
1/+, as depicted in Figures [2^ |2d| Perryman and Wiec- 
zorekl^ observed that the critical value ec equals 4/3 and 
that the connecting orbit is the line 

a; = -^-l (6) 

in the phase plane. For e > W^{S-) and 

change their arrangement, as displayed by Figures 
The unstable manifold W'^{S-) no longer converges to 
the stable node 5'+ such that trajectories from all initial 
conditions close to S- with A > 0 diverge {x{t) Too). 
In this case the parameter A is shifted at a rate that is 
too large for the unstable manifold W'^{S-) to track the 
quasi-steady state Wq. For e > €c but close to €c this 
escape does not occur until A is close to 3 such that one 
would observe the escape only when A is coming to rest 
again, and so there appears to be a lag in the timing of 
escape. 

In summary, for this prototype model rate-induced tip¬ 
ping corresponds to a global bifurcation at parameter 
value ec, a heteroclinic connection from the stable equi¬ 
librium before the parameter ramp to the unstable equi¬ 
librium after the ramp. 


III. DELAY OF EARLY-WARNING INDICATORS AND 
DELAY OF TIPPING 

For the remainder of this paper we will consider the 
following scenario: the speed of the parameter ramp e 
is less than its critical value €c = 4/3 such that without 
noise the system will not tip. The influence of noise, 
which we add to the dynamics of X, will cause the 
system to tip with a certain probability. We can control 


this probability by varying noise intensity and e. We 
choose our parameters such that an escape of x from Wo 
beyond Wq to +oo is extremely unlikely for t far away 
from 0 (and, thus, A far away from Amax/2). We expect 
this escape probability to increase during the ramp of A 
(for t ^ 0). 

The realizations of x for the prototype system 0 
(1^ are governed by the stochastic differential equation 
^DE); 

dXt = [{Xt + \{t)f - l]dt + (7) 


where Wt is standard Brownian motion. The intensity of 
the noise is given by \/2D where D is the diffusion coeffi¬ 
cient. The probability density P{x,t) of the random vari¬ 
able Xt in the SDE Q is governed by the Fokker-Planck 
equation; a linear partial differential equation (PDF): 


dP{x,t) _ d‘^P{x,t) 

dt dx^ 


dx 


f{x,t)P{x,t) 


( 8 ) 


which includes the diffusion coefficient D and drift term 
/(x,t) = {x — 1. Applying Dirichlet boundary 

conditions at some [^Tstart? ^end] will cause the probabil¬ 
ity density to decay over time as realizations escape the 
domain. The probability of escape, Pesc(^n) time step 
tn is therefore defined as: 


Pesc {t n) 


^ _ J P{x,t„)dx 
f P(x,t„-i)dx 


In addition, we use ^ to compute two characteristic 
quantities of the density P{x^t) for 0, the lag-1 auto¬ 
correlation and the variance, shown in Figure These 
quantities are commonly monitored in time series where 
one suspects an underlying parameter drift that ap¬ 
proaches a bifurcation-induced (specifically saddle-node 
induced) tipping point. Both, autocorrelation and vari¬ 
ance, should increase along the time series as the parame¬ 
ter comes closer to its saddle-node value (see Williamson 
and Lentonl^ for other cases such as Hopf bifurcation). 
But what happens for the rate-induced tipping model? 

The lag-1 autocorrelation is defined to be the corre¬ 
lation between successive separated by a time step 

At (we choose At = 0.01): 


_ Cov(X('^_]^;)At5 X-^At) 

""" ^ VVar(X(,_i)At)Var(X,At) 

where X^At is the solution of 0 with density P(-,nAt) 
at time step nAt. 

The initial condition for ^ is the stationary density 
of Q with Ao = A (to) restricted to the fixed domain 
X G^start? ^end] = [~b,2], which Corresponds to the as¬ 
sumption that the ramp-up of A starts from a stationary 
state. See Appendix A(a) for a study of dependence on 
Xend- For this stationary starting point the system can be 
approximately modeled by the Ornstein-Uhlenbeck pro¬ 
cess 
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(a) Autocorrelation 


(b) Variance 


FIG. 3: Traditional early-warning indicators; lag-1 
autocorrelation and variance show a delayed response 
for rate-induced tipping if tipping point is assumed to 
be at the time t = 0 (green dashed line), the closest 
encounter of the stable and unstable manifolds, 
and However, both indicators 

increase before the tipping point if assumed to be at 
t ^ l.b (red dashed line) calculated from the peak of the 
escape rate in Figure]^ Parameters: e = 1.25 and 
D = 0.008 (At = 0.01 for panel 3a) 


corresponds to the closest encounter of the stable and 
unstable manifolds, and lF^(t/+), due to the 

symmetry in the deterministic part. 

Figure displays the lag-1 autocorrelation and vari¬ 
ance for the time interval of most interest, namely t G 
[—3, 2] when system (§ is non-stationary. We ob¬ 
serve that there is a delay in the warning for approach¬ 
ing the tipping point, if we take the tipping point as the 
time t = 0 (green dashed line). The autocorrelation (Fig¬ 
ure has only just started to increase at t = 0. The 
variance (Figure |3b| ) shows an even longer delay in the 
signal. It increases noticeably only after t = 0. Ditlevsen 
and Johnsenll^ concluded for saddle-node induced tip¬ 
ping that only the presence of both indicators, increase 
of autocorrelation and variance, is sufficient evidence for 
the approach of a tipping point. When applied to rate- 
induced tipping, one would conclude initially that the 
warning will be significantly delayed from when we would 
expect the tipping. This warrants a systematic investiga¬ 
tion to see when escape is most likely in close encounters 
with rate-induced tipping. 


dXt = -eXtdt + V^dWt 

where 0 = —/'(—1,0) = 2 is the decay rate at S-. The 
Ornstein-Uhlenbeck process has autocorrelation and vari¬ 
ance given b}i^ 


Autocorrelation: a = exp(—6>At) « 1 — OAt 

D 

Variance: ^ ~ "V 

0 

where we set At = 0.01, thus, giving a = 0.98 and V = 
0.004 in Figure 

We highlight that for starting at to = — oo the sys¬ 
tem will tip with probability one before the ramping 
shift begins. The time of tipping for a stationary sys¬ 
tem is approximated by Kramers’ time, {U{X) = 

-J f(x)dx): 


tk = C exp 


/ AU 


where AU is the height of the potential barrier and the 
prefactor C depends on the curvature at the minimum 
and maximum of the potential. In our case for equation 
0, we have constant values for AU = 4/3 and C = ir. 
We consider the regime where the probability of escape 
from the well, Pesc(^): increases by an order of magnitude 
during the ramping of the parameter A: 

— « maxpesc(t). 

T X 

In this regime we expect the time t for which the escape 
rate Pesc{t) is maximal to occur at the time t = 0. This 


A. Escape rate over time 


To investigate when the escape is likely to occur we 
initially consider the escape rate per unit time calculated 
using the Fokker-Planck equation The escape rate 
per unit time is defined as the fraction of realizations 
that cross a known threshold curve x{t) divided by the 
time step At. We choose a threshold curve x{t) (bright 
blue in Figure]^ beyond which we classify a realization 
as having escaped. The threshold curve is chosen such 
as x{t) = x'^{t) + where x'^U) is the unique trajectory 
of the deterministic part of Q that starts at x(—10) = 
xo = —1 (thus, (x^(t),A(t)) is close to and 

^ = 1.5 is a fixed sufficiently large deviation from x'^{t). 
Appendix |A 2| studies systematically how the choice of 
threshold x{t) affects the results. 


Figure |4| displays the time profile (4a) and the phase 


portrait pb[ ) of the deterministic trajectory x^(t), the 
threshold curve x{t) and the escape rate (over time and 
versus A) obtained via 

In this example, we have chosen e = 1.25, which is close 
to €c = 4/3, and a small noise level D = 0.008. For this 
choice of ramping speed parameter e, tipping would not 
occur in the deterministic case {D = 0). However, with 
a noise level of D = 0.008, roughly 36% of realizations 
that start with initial condition x(—10) = xq go on to 
escape. As in Figure the dark blue curve is the unsta¬ 
ble manifold (the deterministic solution x^(t), 

starting at x(—10) = xq, is extremely close to it). The 
bright blue curve is our threshold curve x{t) (see Figure 
in Appendix |A 2| how the escape time depends on the 
threshold). 

Figure shows that the escape is most likely to occur 
at about t = 1.5, hence, it is delayed, too. The red 
dashed line in Figure [^represents the most likely time of 
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(a) Time profile 
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(b) Phase plane 


FIG. 4: Time profile and phase plane and escape rate 
obtained via for e = 1.25 and D = 0.008. Dashed 
curves represent stable Wq and unstable Wq equilibria, 
the dark blue curve is the unstable manifold 
the bright blue curve is the threshold curve x{t) for 

y = 1-5. 


tipping, given by the peak of the escape rate in Figure]^ 
Therefore taking this as our tipping point we see that 
there is an increase in both the autocorrelation and the 
variance on the approach to tipping. We conclude that 
for this example at least that the tipping is delayed and 
thus the early-warning signals can still offer forewarning 
of rate-induced tipping. 

For a systematic study of how the time of most likely 
escape depends on the system parameters we formulate 
a variational optimization problem for the optimal path 
of escape. 


IV. MOST LIKELY (OPTIMAL) ESCAPE PATHS — 
THE GENERAL VARIATIONAL PROBLEM 


In this section we will formulate the ODE boundary- 
value problem (BVP) determining locally most likely 
paths for escape. 

We define the most likely escape path as the path going 
from given xq to a given xt in a time interval [to, Tend] 
maximizing the functional 


F = exp 


Uo-Ut 

2D 


r 

Jto 


Tind / ^2 


AD 


+ V ]dt 


along the path. The terms in F are 
U{x,t) = - j f{x,t)dx, 


Vs{x,t) = 


- J_t^V _ _ J_^ 


AD \ dx 


2 dx‘^ 2D dt ’ 


( 10 ) 


( 11 ) 


For a differentiable path x the functional F equals the 
probability of a realization Xt following a sequence of in¬ 
finitesimally small intervals [x{kAt) — (^/2, x{kAt) F S/2] 
in the limit 0 < S At 1 (up to a constant factor in¬ 
dependent of x). Recall, that the random variable Xt had 
a probability density function P(x, t) given by the linear 
Fokker-Planck equation ® , from which the functional F 
is derived. See Zhangl^,^. 25-31) for a detailed deriva¬ 


tion of the functional F from equations (10)-(12) for a 
time independent potential U{x) (a simple extension is 
made for a time dependent potential U{xD) in Lin and 
Ho^^ and Ho and Dai^^ . 

Assuming a fixed time interval [to, Tend] and fixed start 
and end points xq and xt, local critical points of F 
are given by the Euler-Lagrange equation, a 2nd order 
BV^ 


X = 


x{tQ) = a;o, 
^(Tend) ~ ^T- 


(13) 


We would like to point out that the BVP ( [l^ used 
to calculate the locally optimal path is valid for a scalar 
time-dependent system and for finite (non-small) noise 
variance 2D. In the small noise limit, one can use min¬ 
imum action methods to find the optimal path, which 
can be applied to multiple dimension^^. Eurthermore, 
according to Ren, Vanden-Eijnden et in gradient 

systems, over an infinite time interval, the optimal path 
becomes a minimum energy (where ‘energy’ refers to the 
functional F that is optimized) path that forms a het¬ 
eroclinic orbit between the two local minima of the po¬ 
tential. However, Eigure demonstrates that even for 
relatively small noise levels D (such as D = 0.008 as cho¬ 
sen for previous illustrations) we are far away from the 
small noise limit such that Tend is of order 1: Eigure [5a| 
shows the (locally) optimal path x{t) for Tend = 20. We 
observe that for a long time (1 < t < 18) the path x{t) 
stays close to the saddle [/+ before eventually escaping 
to the chosen xt = A. As Eigure \Sh\ shows, the lingering 
of x{t) close to the saddle is only optimal for the fixed 
large Tend = 20. The functional M = log(F) increases 
for decreasing Tend- 

This implies that for positive (even small) noise vari¬ 
ance 2D the functional F should also be optimized with 
respect to the traveling time Tend of the path. We formu¬ 
late the extended BVP corresponding to critical points 
with respect to path and traveling time in rescaled time 
on the base interval [0,1]. The BVP will then be solved 
with standard continuation software AUTCj^. The BVP 


(13), rescaled to [0,1] is (split into two components): 


Vo = V(xo,to), Ut = V(xT,Tend)- 

The quantity U is the potential of the deterministic 
part / of the SDE ([^ such that ([^ can be written in 
terms of U (x, t): 


±1= X2{Tend-to), Xi{0) = Xq, (14) 

X 2 = 9{xi,t){Tend - to), Xi{l) = Xt (15) 

where to (fixed) and Tend (free) are the start and end t 
values and 


dXt = + V^dWf (12) 

OJit 


g{xi,t) = 2D 


dx 


{xi{t),t). 
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The solution of (14)-(15) is a critical point of F, given in 
( 10 ), among all possible paths connecting from xq to xt 
in a fixed time T = Tend ~ ^o- The function M = log(T), 
written for the rescaled path is 


M = 


i 

/ 


Uo-Ut 

2D 


AD 


Vs{xi{t),t) (Tend - to)dt 
(16) 

where UQ = U{xo,to), Ut = U{xT,Tend)- 
Paths maximizing F (and, hence, M) also maximize 



(a) (b) 


Therefore, have to solve the four-dimensional BVP (14), 
0 for xi(t), X 2 {t)_,zi{t)^ with the addi¬ 

tional integral condition (18) and the additional param- 
eter Tend- We use AUTO (Version: AUTO-07P)P to 
study the solutions of 0 , 0 , 0 in depen¬ 

dence of the system parameters D and e. 


V. SEQUENCE OF CONTINUATION STEPS FOR THE 
OPTIMAL PATH TO ESCAPE IN OPTIMAL TIME 


Since (14), (15), (17), (18) is nonlinear we need a se¬ 
quence of initialization steps to arrive at the optimal path 
for particular desired values of ramping speed parameter 
(initially e = 1.25, close to critical value ec = 4/3) and 
noise variance (initially 2D = 0.1). An advantage of us¬ 
ing continuation is that once we have obtained an optimal 
path in an optimal time for a particular set of parame¬ 
ters we are free to perform a systematic parameter study 
of solutions of (14), 0 , 0,0 varying the ramping 
speed e and noise level D. 


FIG. 5: (a) Optimal path for Tend = 20. (b) Plot of 
function M that needs to be maximized w.r.t. Tend- 
to = -10, e = 1.25, D = 0.008. 


the probability of realizations of SDE ([^ following it. 
Figure 5b plots M along paths satisfying (p^-(p!^ for 
a range of end times Tend- Us maximum corresponds 
to the time Tend for which the functional M is (locally) 
maximal among the range of Tend shown (this is for the 
fixed positive but small noise variance 2D = 0.016). We 
now extend the BVP (14)-(15) to include the criticality 
of Tend into the optimization problem. We outline the 
BVP of the variational problem for the general case here. 
The specific case for this example is in Appendix [B| 
Introducing the derivatives of xi and X 2 w.r.t. Tend as 


Zi{t) = 


dxi{t) 


aTend ’ 
these derivatives satisfy 
Z\ — X2 T -2)2 (Tend fo): 


Z2{t) = 


dX2{t) 

5Tend ’ 


u(o) = o, 


^2 = 5(a;i,d +-^|^/^u(Tend - ^o), 2 : 1 ( 1 ) = 0^ ^ 

Critical points of M(xo, Iq^xt^ Tend 5 xi{'),X 2 {')), given in 
equation (16), w.r.t. Tend satisfy 

dM dM dxi dM 8x2 

which produces the integral condition: 


0 = m := 




1 dU{xT,Tend) , X2{tf 


2D STend 


AD 


^ f x2it)Z2{t)+g{Xi{t),t)Zi{t) d-^o) 


Vs{xi{t),t) 

dt (18) 


A. List of free parameters 

First we discuss some of the parameters used and rea¬ 
soning for their initial values, as given in Table [ij We 

TABLE I: Types of parameters used in continuation 
steps and their initial values 


System (fixed) Continuation 

Bifurcation 

Monitoring 

parameters 

parameters 

parameters 

parameter 

p= 1 

Tinit = 0 

1.25 

M = Mo 

Amax — 3 

XT = Xo 

D = 0.05 


to = -10 

Tend = — 9 



Xo = -1 

m = mo 




introduce the factor Tinit in equation (15) as an artificial 
parameter. Thus, 0-0 changes to 


Xi — (Tend fo) 

^2 ^(^l:^)(Tend fo)Tinit5 


giving a trivial system {xi = X 2 (Tend “ ^2 = 0 ) for 

Tinit = 0 ? connecting it to system ([T4|)-(15) via a con¬ 


tinuation in Tinit from 0 to 1. Eurthermore, we initially 
choose Xt = xq and Tend = fo + 1 = close to the 
initial time value to. Einally, the parameters M and m 
are used to monitor the values of the integrals in equa¬ 
tions ( [T^ and ( p^ , respectively, such that sign changes 
of m correspond to critical values of the functional M 
(and, hence, T). The initial values of M and m are set 
equal to the integrals in (H!^ and (18) along the initial 


path. Outlined in TablehllTs a brief summary of the con¬ 
tinuation steps performed to create an optimal path in 
an optimal time. We proceed with a brief explanation 













































for each of the continuation steps with a more in depth 
discussion provided in Appendix [C] 


TABLE II: Summary of continuation steps to perform 
in order to achieve an optimal path for escape in an 
optimal time (in brackets are the values used) 


Step 

# 

Continuation 

parameter 

Initial value 

End value 

Other free 
parameters 

1 

Tinit 

(0) 

(1) 

m, M 

2 

Xt 

Xo (-1) 

> 1 (4) 

m, M 

3 

Tend 

~ to (—9) 

» 1 (20) 

m, M 



B. Step 1: continuation of Tinit from 0 to 1 


The first step ends in a solution of the full system of 
equations (14)-(15). This is the short orbit (blue) from 
xi = S- = —1 to xi = —1 in Figure The param¬ 
eters M and m are kept free during the continuation, 
monitoring the integrals in (16) and (0^81). 


FIG. 6: Illustration of trajectories (colored) at different 
stages of the xt continuation step superimposed on the 
phase portrait for the full rate-induced system, when 
stationary, A = 0. xt = xq = —1 (blue), xt = 0 (red), 
Xt = 1 (green), xt = 4 (black). 


VI. OPTIMAL PATH FOR ESCAPE FOR NOISE AND 
RATE-INDUCED TIPPING 


C. Step 2: continuation in xt 


The result of Step 1 is a path maximizing M when trav¬ 
eling from xo = — 1 to Xt = — 1 in unit time (to = —10, 
^end = —9), where A(t) is close to stationary (A ^ 0). 
Step 2 changes the right boundary value xt to the de¬ 
sired location. Guided by our aim to find most likely 
paths for escape, we perform a continuation to xt = 4 
(we count a trajectory of 0 that reaches xt = 4 as 
having escaped). Figure]^ shows a sequence of solutions 
of ([l4|)-(p!^ (colored) for the continuation stages of xt 
superimposed onto the phase portrait of the system (14)- 


(15) for A = 0. The results of this step is a path maximiz¬ 


ing M that connects xq = — 1 and xt = 4 in a very short 
time period (Tend — to = 1, see Figure 14b in Appendix 
G2). Glearly this is not the optimal time to make this 


transition and so the next step is to continue in Tend to 
get to more realistic timings of escape. 


D. Step 3: continuation in Tend 

We increase Tend to a large value, monitoring M and 
m. Figure |5l^ shows the graph of M over Tend, which has 
a pronounced maximum at Tend ~ 1-5. Gritical points 
of M are detected when m changes sign, and therefore 
gives the optimal path for an optimal time provided M 
is a maximum. 

The optimal path constructed by the above steps will 
be systematically continued in the system parameters D 
(noise variance 2D) and e (ramping speed of A(t)) in Sec¬ 
tion IVlIl 


This section compares the optimal path for escape with 
the escape calculated directly from the solutions of 0. 
We are interested in the timing of escape, defined as the 
timing of crossing certain threshold curves. We do not 
want to have to rely on running full Monte Garlo simula¬ 
tions, but instead use the optimal path theory developed 
by Zhangl^, Ghaichian and Demichevl^. We include 
the optimal path for escape (in green) into Figure]^ to 
compare the most likely timing of escape (the peak in 
the escape rate, measured at the threshold x) with the 
time t the threshold x intersects with the optimal path as 
computed through continuation, see Figure The cor¬ 
responding phase portrait in Figure [7b| shows that both 
the simulations and optimal path suggest that the escape 
does not happen until just short of A = 3, the moment 
the potential well or steady states Wd, Wq are coming 
to a rest. 

Figure [^illustrates that the optimal path matches the 
mode (peak) of the escape rate well. In general, if the 
escape rate over time is unimodal with a sharp peak then 
the time profile of the optimal path is a good description 
of this peak. More precisely, the mode of the escape rate 
occurs very close to a time t for which x = Xi(t), where 
xi is the first component of the optimal path, the solution 
of the extended BVP (14), ( [T^ , ( [T7| ), ( [T8| ). 


A. Dependence on choice of threshold curve x{t): 

The choice of the threshold curve x(t) in Figure ^is at 
x(t) = x^(t) with ^ = 1.5 (recall that x^(t) is the tra¬ 
jectory of the deterministic part of (j^. Figure]^ shows a 
color plot for the escape rate at the threshold depending 
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t 





(a) Time profile. 


(b) Phase plane. 


FIG. 7: Time profile and phase plane with escape rate 
from simulations (top) and optimal path added in green 
for e = 1.25 and D = 0.008. Dashed curves represent 
stable Wq and unstable Wq quasi-equilibria, the dark 
blue curve is the unstable manifold the bright 

blue curve is the threshold curve x{t) = + y with 

y = 1-5. 


on the distance y. The distance of the optimal path from 
the unique trajectory x'^{t) added in white. This high¬ 
lights that, provided the threshold x{t) is sufficiently far 
from x^(t), the optimal path will cross the threshold at 
the same moment the escape rate through the threshold 
is at its peak. 



t 


FIG. 8: Golor plot of the escape rate at the threshold 
x{t) depending on the distance y from the unstable 
manifold W'^{S-). Distance of optimal path from x'^{t) 
added in white, e = 1.25, D = 0.008. 


B. Dependence on starting time to and starting position 

xq: 


We emphasize that the optimal path calculates the lo¬ 
cal optimum, i.e. assuming the system has not tipped be¬ 
fore the ramping shift begins. We set to = —10, xq = — 1 
to represent starting at the bottom of the potential well 
at time —oo. Ghanging the starting time, for example 
to to = —15 or to = —5 has no effect on when the opti¬ 
mal path escapes and only extends or shortens the time 
profile of the path presented in Figure [^a). Likewise, 
changing the starting position, provided it is still inside 
the well, has no effect on when the optimal path escapes. 
For slightly different xq the optimal path will converge 
onto the path in Figure [^a) before t = 0 and then follow 


the same path for escape. 

Figureshows how the xi{t) component of the opti¬ 
mal path is connected to the evolution of the density of 
realizations. 


High 



FIG. 9: Time profile and phase plane for density plot of 
simulations and optimal path added in bright blue for 
e = 1.25 and D = 0.008. Dashed white curves represent 
stable and unstable equilibria. 


In the time profile plot, initially the spread of the dis¬ 
tribution is very narrow centered around the steady state 
Wo, due to the small noise level D. When the system 
shifts, this distribution widens reflected by a lower den¬ 
sity over a larger x range. Once the shift stops, the den¬ 
sity gradually becomes concentrated again, but some re¬ 
alizations have escaped, indicated by the elevated density 
at X = 4. Initially, the optimal path is right at the mean 
of the distribution. The time when the optimal path de¬ 
viates from the mean equals the time the density in the 
simulation is at its widest and where the additional mode 
(at X = A) appears. This once again suggests that the 
optimal path, derived from BVP (14), (15), (17), (18), 
describes the escape of realizations of the stochastic dif¬ 
ferential equation 0. 


VII. TIMING OF ESCAPE IN 2 PARAMETER PLANE 


One of the advantages of reducing the study of escape 
time to optimal paths is that we can perform a system¬ 
atic parameter study with moderate computational ef¬ 
fort. First, we investigate the timing of escape in the two 
parameter plane of the ramping speed e and noise level D, 
Figure panels |lQa[ |lQc| and |lQe[ Panel |lQa| indicates 
with a black marker the end time Tend at which a partic¬ 
ular optimal path reaches the end position xt = A. Panel 
|lQc| shows a color contour plot, with the color denoting 
the time Tend foi* a range of optimal paths dependent 
on the ramping speed e and noise level D. Recall that 
the optimal path is calculated by solving the system of 
equations (14)-(18) and following the continuation steps 
outlined in Section |V| for a particular e and large noise 
level D. Then for each e a final continuation is performed 
over D to create an 11 x 40 grid for the color plot. 

The optimal path begins at to = —10 and so the length 
of the time interval for the path is between 11.7 and 13.5 
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time units for this range of e and D values. This demon¬ 
strates that for a small e value and small noise levels 
escape occurs for positive Tend? that is, with a delay com¬ 
pared to the time of closest encounter of the determinis¬ 
tic manifolds and (which would be at 

t = 0). As e increases towards Cc = 4/3 and the noise 
level increases the time to escape decreases. This can be 
seen more clearly in panel |lQe[ which takes cross sections 
of panel |lQc| for different values of e. The relationship 
is almost linear between the logarithm of the noise level 
and the time at which the final destination xt is reached. 
In summary, panels and indicate that the escape 
occurs with a delay especially for small noise. To inves¬ 
tigate the precise value of the delay we will look at the 
timing of intersection between the optimal path and the 
stable manifold indicated by the black marker 


in panel 10b 


The reason for considering the stable manifold 
as a threshold is that it plays a role similar to a saddle 
in stationary escape problems. Once a realization has 
crossed this threshold it is more likely to escape to Too 
(in finite time in our example). One may expect the op¬ 
timal (that is, most likely) escape path to cross this man¬ 
ifold when the two manifolds W'^{S-) and are 

closest together at t = 0. The question is then whether 
the escape across the stable manifold occurs at 

t close to 0. Panels |lQd| and |lQf| present the timing of 
crossing the stable manifold to establish if this 

is the case. 

We observe that the range of crossing times is smaller 
than the range of end times to reach xt{= 4). This is 
expected since the traveling time from to xt de¬ 

creases for increasing noise level D. For small noise levels 
the optimal path tracks the manifold for longer. 

In the limit of large noise level D in the (e, T))-parameter 
plane the most likely crossing time is tcross ~ to (not 
shown). In this limit we have a purely noise-induced tran¬ 
sition as the potential is nearly stationary at to- For de¬ 
creasing D the intersection between the optimal path and 
stable manifold IT^(t/+) varies with different ramping 
speeds e such that we have combination of noise and rate- 
induced tipping, with timing depending on both parame¬ 
ters. As the ramping speed e and noise level D decreases 
the crossing time delay tcross increases. For the smaller 
noise levels in the computed range the intersection tcross 
of the optimal path with the stable manifold IF^(I/+) is 
of order 1, when the manifolds and IF^(t/+) are 

significantly further apart than at t = 0. This justifies 
the claim in the abstract that for noise- and rate-induced 
tipping the escape is delayed in the small noise limit. 
Figure gives a crude estimate for the probability of 
escape depending on the ramping speed e and noise level 
D. Note though, these are not the true probabilities but 
rather the values of the functional M = logT (see @) 
which the most likely path optimizes. The color in Fig¬ 
ure 11a in the 2 parameter (e, T))-plane equals the value 
of the functional M along the optimal path found at the 
corresponding point in the (e, T))-plane. As expected, the 



FIG. 10: Plots for end time Tend (^)?(c),(e) and the 
crossing time tcross with the stable manifold IF^(I/+) 
(b),(d),(f) of the optimal path (solution of (p^-(p^) for 
escape. (a),(b) Time profile of optimal path (green) and 
stable manifold (red) with the black marker 

highlighting Tend (a), Across (b) for a particular optimal 
path (e = 1.25, D = 0.008). (c),(d) Color contour plots 
for end time Tend (c) and the crossing time tcross (d) in 
the 2 parameter (e, T^) - plane. (e),(f) Cross sections of 
(c),(d) respectively where each contour represents 
different value of e, spaced evenly at 0.04 intervals, 
starting with e = 1.05 (dark blue, top) increasing to 
e = 1.25 (bright blue, bottom). 


largest probability of escape is for large ramping speeds 
and large noise levels. The value of M is smallest for 
slow ramping speeds and low noise levels. Figure |llb| 
displaying the cross section of Figure |lla| for different 
values of e illustrates that M decreases logarithmically 
as D decreases on a logarithmic scale. 


VIII. GENERAL DELAY OF TIPPING 


In Section [Vnj we have shown that the tipping is de¬ 
layed especially for small noise levels. In the context of 
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FIG. 11: (a) Color plot of the value of M. (b) Cross 
section of (a) each contour represents different value of 
e, spaced evenly at 0.04 intervals, starting with e = 1.05 
(dark blue, bottom) increasing to e = 1.25 (bright blue, 

top). 


autonomous systems Bakhtin^^ gives an asymptotic for¬ 
mula for this delay. Bakhtinl^ considers rare escapes for 
small noise levels for a process dx = xf{x) -f \/ 2 DdWt 
on the interval [^4, B] containing 0, starting from xq < 0, 
where / is uniformly positive. Then the first time 
T = Tend — ^0 to exit at the point B (under the condition 
that x{t) does indeed exit at B) in the limit \/2T ^ 0 
satisfies 


T = Cl In 


VW 


where ci is a constant independent of D. This states 
that for an autonomous system, the time for rare es¬ 
capes increases linearly as the noise level decreases ex¬ 
ponentially. This is consistent with our findings for the 
non-autonomous system, that as the noise is decreased 
the time tcross ^.t which the optimal path crosse s the sta¬ 
ble manifold W^{U-^) increases slowly. Figure [TQf[ To 
conclude, we find a similar relationship for the delay in 
the rate-induced tipping as that of Bakhtin^^ for rare 
escapes of an autonomous system. The observed level 


of delay in Section VII is of order 1 such that the noise 
levels that we consider small in Section EU are still far 
larger than the small-noise limit, for which BVPs for op¬ 
timal escape paths are available in arbitrary dimensions. 
(These paths tend to be connecting orbits such that the 
optimal time is always infinity)!^. 


We extended the boundary-value problem for the most 
likely path for tipping (escape) based on Zhangl^ to in¬ 
clude optimality of time for finite noise. This additional 
optimality criterion created a variational optimization 
problem that we solved computationally with continu¬ 
ation techniques (using the package AUTO). With the 
help of continuation we performed a systematic param¬ 
eter study in the (e, T)-plane (ramping speed vs. noise 
level). The time when the optimal path for escape crosses 
the stable manifold 1U^(1/+) is a measure for the timing 
of tipping. We find that for large ramping speeds and 
noise levels there is no delay and even for lower ramping 
speeds there is only a small delay. However, for small 
noise levels D the tipping delay is of order 1. 

We hypothesize that the observed delay in tipping is 
present independent of the particular form of X{t) as long 
as it is qualitatively similar to the ramp like shift (§. 
Similarly, this delay should be observable independent of 
the particular shape of the potential well l/(',t). This 
paper demonstrated that the optimal path for escape, a 
solution of a BVP, matches simulation results well. The 
technique used to find the optimal path of escape finds 
the local maximum and is general such that it can be 
used to determine the timing for any type of tipping. 

However, the optimal path may miss the global opti¬ 
mum when there is more than one realistic opportunity 
for escape. For a small single window of escape as consid¬ 
ered in this paper the escape rate will form a unimodal 
distribution with a narrow peak, for which the optimal 
path is close to the mode. However, if one considers dif¬ 
ferent scenarios X{t) for ramping the system parameter 
(for example, one that is not monotonically increasing, 
see Ashwin, Perryman, and Wieczorek®) , the escape rate 
would have a multimodal distribution. We conjecture 
that we find one optimal path for each of the modes of 
the distribution. It is unclear if for non-monotone param¬ 
eter shift A(t) the tipping or the early-warning indicators 
are delayed for the small noise levels. This would fur¬ 
ther support the conclusion that the autocorrelation and 
variance can be used as early-warning signals for rate- 
induced tipping events. Furthermore, this paper has fo¬ 
cused on the one-dimensional case. Thus, an extension 
to the general multiple dimensional case is still required. 


Appendix A: Dependence on parameters 


IX. CONCLUSIONS 

We have shown that two commonly used early-warning 
indicators of tipping (increase of autocorrelation and in¬ 
crease of variance) are present but delayed in a proto¬ 
typical model of rate-induced tipping. By looking at the 
timing of escape using optimal paths we find that the 
tipping event itself is delayed for small noise levels. We 
conclude that the delay in the early-warning indicators is 
consistent with the delay in the actual tipping (at least 
for the example). 


This appendix details how the choice of the parame¬ 
ters Xend: the Upper domain boundary, and the threshold 
parameter y affects the results presented in the paper. 

1. Domain boundary parameter 

We investigate the effect the upper boundary of the 
domain, Xgnd, has on the early-warning indicators, the 
increase of autocorrelation and variance. We choose a 
domain [xgtart, ^end] with Xgtart = -6 and Xend as shown 
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in Figure Choosing the domain fixed in time is nat¬ 
ural as in realistic problems we do not know a-priori the 
location of the moving well. In Section III we chose 
^end = 2, which corrcspouds to a wide domain for the 
problem (including realizations into the computation of 
the early-warning indicator, which are already escaping). 
Figure explores the effect narrowing the domain has 
on the decay rate estimate 0 = (1 — a)/At (where a is the 
lag-1 autocorrelation with time step At) and variance V. 
We consider the decay rate, instead of the linearly related 
autocorrelation as the decay rate is independent of the 
time step At. 



t t 


(a) Decay rate (b) Variance 

FIG. 12: Effect of the width of domain on the decay 
rate and variance by varying upper boundary Xend 


In Figure |12a| we see that the timing of the onset of 
the decrease of the decay rate estimate 0 is independent 
of the upper boundary of the domain, Xgnd, and hence, 
so is timing of the onset of the increase of the lag-1 au¬ 
tocorrelation a. Likewise, in Figure |12b| the timing of 
the onset of the variance is independent of Xend and, im¬ 
portantly, shows no increase before t = 0. The precise 
values of the autocorrelation and the variance depend, of 
course, strongly on the width of the domain (and, hence, 
on Xend)- Thus, wc have shown that, while autocorrela¬ 
tion and variance change quantitatively with the domain 
width, the timing of their increase (which is the early- 
warning indicator) does not change. 


2. Threshold parameter 



t t 


(a) D = 0.1 (b) D = 0.008 

FIG. 13: Evaluating distance y required between the 
deterministic trajectory and threshold x{t) in 

both the large and small noise limit cases, e = 1.25 


critical value ~ 1 ? large fraction of realizations would 
count as escaped even for t close to ±10, where A is close 
to stationary, in the large noise case. Eor smaller noise 
level Eigure |13b| illustrates that a smaller y will suffice. 
Both figures show that above a certain minimal value of 
y the window of escape remains nearly independent of y. 
The value y = 1.5, used in Section [Vl] for our threshold 
curve x(t), is well above that minimal value of y. 


Appendix B: Variational optimization problem for specific 
example 


The following is the variational optimization problem 
for the specific rate-induced example discussed in the pa¬ 
per, equation 0 - 

The potential U{x^X{t)) for equation 0 is: 

^3 

U{x, A(t))= — —— Xx^ + (1 ~ such that 

o 

U'= -a;2 - 2Aa; + 1 - 
U"= -2{x + A) 

U= — Xx(^x “h 2A) = —€Ax(Aniax — '^) (^ 2A) 


wuere u and u represent wx ^xx^ p^u^xx- 

tial w.r.t. space and time respectively and A is given by 
equation 0 These equations feed into the W, equation 
O giving: 


This section presents in more detail how the distance 
y of the threshold curve x{t) (at which we consider a 
realization as escaped) from the deterministic trajectory 
x'^{t) influences our results {x{t) = x'^{t) ± ^). If we 
choose y too small, then escape will be detected every¬ 
where. This is demonstrated by Eigure Eor example, 
when the steady state is stationary no escape should be 
detected, because, if realizations cross the threshold x{t)^ 
most will not escape to ±oc but will return to the unsta¬ 
ble manifold W'^{S-). Clearly for a larger noise level D 
a greater y is required as there will be larger fluctuations 
about the unstable manifold than for a small 

noise level. 

Eigure [T^ demonstrates that for values of y less than a 


± 4X^x^ ± (1 - Xy ± 4Ax3 - 2x(x ± 2A)(1 - A^) 
" “ TD 

^Ax(Aniax A)(x ± 2A) 

+ -2C- 

and so the 2nd order boundary value problem, equation 
split into 2 first order ODEs augmented with the 
ODE for A, equation 0 , which are to be solved on the 
[0,1] time domain in AUTO looks like: 

±1= X2(Tend - to) (Bl) 

±2= h2(xi,A(t))(Tend “ to) (B2) 

A= /l3(A(i))(Tend — to) (B3) 


Escape rate 
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where (B1)-(B3) correspond with (14)-(15) and 


^ 2 (^ 1 , A(t))= 2D 


dVs{xi,X{t)) 


dxi 

h 3 {X{t))= eA(Amax - A) 


The function M = log(F), equation ( [l^ for the general 
case, is maximized and used to monitor any maxima or 
minima, is given by: 


M = 


/[ 


U{xo, A(to)) - U{xt, A(Tend)) 
2D 


■(^+V(a:i,A(f)))(Tend-to) 


The variational equations for zi , Z 2 0 and Z 3 = 
are given as: 


df (B4) 
d\{t) 


Z\ —X2 + ^2 (lend — to) 

Z2 =h2{xi,X{t)) 

' dh2{xi,X{t)) 


(B5) 

(B6) 


9a; 1 




+ (B7) 

To locate the local maximum of M, which is the deriva¬ 
tive of (B4) w.r.t. Tend we have a second integral condi¬ 


tion corresponding to equation (18) in the paper (multi¬ 
plied by —AD to remove fractionSp 


L 


Z2 + 2h2{xi, X{t))zi + j 


dU{XT,X{Tend)) I 4 mn-r \(+\\ 

+2- ^ -ha;2 +4Z;Vs(a;i, A(f)) 

oiend 


df = 0 (B8) 


and thus for the general example in the paper which had 
five equations to solve (14 )-([T^, 0 and and for the 
specific r ate- i nduc e d ex a mple ther e are seven equations 
to solve (|B1|)-(|B3|), (|B5|)-(|B7|) and 


Appendix C: Detailed explanation of continuation steps 

We provide further explanation of the continuation 
steps pre sented in the paper for the specific rate-induced 
example (|B1|)-(|B3|). 


having started with Tend = —10 and Tinit = 1- Starting 
with Tend = —10 would mean: 


xi = 0 xi{— 10 )= —1 

X 2 = 0 xi{—9) = —1 

which has no unique solution. Incorporating the artificial 
continuation parameter Tinit (initially at 0) and setting 
^end = — 9 we have: 


xi = X 2 xi{— 10 )= —1 
X 2 =0 xi{—9) = —1 

which does have a locally unique solution. Thus, we can 
continue in Tjnit until Tinit = 1 to obtain a solution of the 


dT^nd system (B1)-(B3). Note that the parameters M and 


m have to be kept free d urin g this continuation such that 
the integral conditions (B4), (B8) are always satisfied. 


2. Step 2: xt continuation 

The sketch of the phase portrait. Figure is an ac¬ 
curate representation of the full system that does not 
change over the time considered, since this continuation 
is for to = —10 and Tend = —9 and hence A « 0. The 
phase portrait contains two saddles, which are located 
close to the equilibrium points of S- and U- and one 
center close to the origin. The saddles are offset to the 
left of S- and U- by approximately D/ 2 . The center of 
the elliptic region is located at (xi,X 2 ) ~ (T,0). 

The initial trajectory after step 1 {xt = xq = — 1) is 
a unique solution that is contained within the elliptic re¬ 
gion near the saddle (see blue trajectory Figure]^. The 
boundary condition xi{to) = xq ensures that all trajecto¬ 
ries start on the dashed line S- = —1. During continua¬ 
tion in Xt the trajectories need to travel further but still 
in the same time interval (to = —10, Tend = ~9) as xt 
is increased. Therefore, the starting position increases in 
the X 2 direction where the vector field has a larger xi 
component. This enables the trajectories to travel faster 
in order to travel further in the same time period. See 
Figure for intermediate phase portraits in this contin¬ 
uation, Figure pTa| for the phase portrait and |14b| for the 
time profile of the (final) trajectory with xt = 4. 

This trajectory corresponds to the optimal path for a 
purely noise-induced tipping since A is close to stationary. 
One would expect the optimal time for tipping to be a 
result of both noise and rate-induced tipping. 


1. Step 1: Tinit continuation 

This continuation is similar to an integration in time 
continuation. However, there is a difference between this 
continuation and performing a continuation in Tend to —9 


3. Step 3: Tend continuation 

The next step is to perform a continuation in Tend while 
monitoring m for roots (or M for critical points). Since 
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(a) = 4 



t 


(b) Optimal path after xt 
continuation step. 


FIG. 14: (a) Trajectory in (xi,X 2 ) - plane, (b) optimal 
path after xt continuation. to = —10, e = 1.25, 

D = 0.05. 


M may have several local minima and maxima for in¬ 
creasing Tend we continue in Tend too sufficiently large val¬ 
ues where we observe the asymptotic monotone decrease 
of M. In our example, we continued Tend from —9 to 
20, monitoring the bifurcation diagram in the (Tend,^) 
- plane (not shown but similar to Figure [51^. 

For choice of parameter values displayed in Figure |5l^ 
there is only one critical point of M, which corresponds to 
the maximum we are interested in. Figure |15a| displays 
the optimal path after the Tend continuation step, for 
^end = 20. For Tend = 20 the optimal path reaches the 
saddle at t ~ 1, but waits at the saddle until t ~ 18 
before escaping to xt, which is optimal only in the limit 
D ^ 0. Therefore, we dete ct when m = 0 which satisfies 
the integral condition (B8) to identify the maximum of 
M and hence achieve the optimal path in an optimal 
time, see Figure [T^ 



t 

(a) Tend = 20. 



t 

(b) Tend ~ 1.43. 


FIG. 15: Gomparisen between optimal paths for 
^end = 20 and after m continuation is completed, 
to = -10, e = 1.25, D = 0.05 
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